#Replication material for Goet, Niels D., "The Politics of Procedural Change."
#First version: 1st May 2018
#This version: 17th February 2019

# dependencies (please uncomment the lines below to install the required R dependencies)
# dependencies <- c('ggplot2',
# 'bcp',
# 'strucchange',
# 'foreign',
# 'optimx',
# 'lme4',
# 'grid',
# 'stargazer',
# 'spduration',
# 'survival')

# lapply(dependencies,function(x) install.packages(x))

######
# Mixed generalised linear model (binomial)
######
scriptPath <- rstudioapi::getSourceEditorContext()$path
scriptPath <- gsub("runReplication.R","",scriptPath)
setwd(scriptPath)

UKHCSOdat <- read.table("data/UKHCSOdat.tab",sep='\t',header=T)

UKHCSOdat$daysUntilElect <- as.Date(UKHCSOdat$dateNextElection) - as.Date(UKHCSOdat$state_opening)

require(lme4)
library(ggplot2)

####################
# Pre-processing
####################
workloadVars <- subset(UKHCSOdat,select=c(session_ref,totalMotions,parl_size,session_days))
prin_comp <- prcomp(subset(workloadVars[complete.cases(workloadVars),],select=-session_ref), scale. = T,center=T)

workloadDat <- data.frame(workload=data.frame(prin_comp$x[,1])[,1],session_ref = workloadVars[complete.cases(workloadVars),]$session_ref)

UKHCSOdat <- merge(UKHCSOdat,workloadDat,by="session_ref",all.x=T)

UKHCSOdat$daysUntilElect <- as.numeric(UKHCSOdat$daysUntilElect)

# Rescale variables
####################
reScale <- function(var){
  var <- (var - mean(var,na.rm=T))/(2*sd(var,na.rm=T))
}
UKHCSOdat$maj_size <- reScale(UKHCSOdat$maj_size)
UKHCSOdat$obstruction <- reScale(UKHCSOdat$obstruction)
UKHCSOdat$daysUntilElect <- reScale(UKHCSOdat$daysUntilElect)
UKHCSOdat$sdGov <- reScale(UKHCSOdat$sdGov)
UKHCSOdat$sdOpp <- reScale(UKHCSOdat$sdOpp)
UKHCSOdat$govOppDiff <- reScale(UKHCSOdat$govOppDiff)
UKHCSOdat$partyStrength <- reScale(UKHCSOdat$partyStrength)
UKHCSOdat$workload <- reScale(UKHCSOdat$workload)


# Run GLM models                    
library(optimx)
library(lme4)
mod1 <- glmer(dv_bn ~ 
                obstruction +
               sdGov +
               govOppDiff +
                
                workload +
                change_control + 
                partyStrength +
                maj_size + 
                daysUntilElect +
                reformprevsession + 
               (1|reformID), data=UKHCSOdat, family=binomial,
control = glmerControl(optimizer = "bobyqa"),
nAGQ = 10) 

mod2 <- glmer(dv_bn ~ 
                obstruction +
                sdOpp +
                workload +
                change_control + 
                partyStrength +
                maj_size + 
                daysUntilElect +
                reformprevsession +
                (1|reformID), data=UKHCSOdat, family=binomial,
              control = glmerControl(optimizer = "bobyqa"),
              nAGQ = 10) 


# Generate coefficient plots
statSigCol <- "black"
statInsignCol <- "grey80"
ratioFormat <- T
xlab <- "Log odds"

model <- "GLMM"
source("customFunctions/genCoefplot.R")
mod1Plot <- genCoefplot(mod = mod1,
            leftAxis = T,
            statSigCol = statSigCol,
            statInsignCol = statInsignCol,
            plotTitle = "Model 1",
            xlab = xlab,
            model = model)

mod2Plot <- genCoefplot(mod = mod2,
                        leftAxis = F,
                        statSigCol = statSigCol,
                        statInsignCol = statInsignCol,
                        plotTitle = "Model 2",
                        xlab = xlab,
                        model = model)

#generate coefplot
##################
library(grid)

pdf("figures/FIG1_GLMM_COEFPLOT.pdf",width =.8*8, height = .8*3)
grid.newpage()
grid.draw(cbind(ggplotGrob(mod1Plot), ggplotGrob(mod2Plot), size = "last"))
dev.off()

#write out model file
library(stargazer)
mod_stargazer <- function(output.file, ...) {
  output <- capture.output(stargazer(...))
  cat(paste(output, collapse = "\n"), "\n", file=output.file, append=F)
}
mod_stargazer("models/GLMMModels.tex",mod1,mod2,
          model.names=F,
          label="GLMMmod",
          title= "GLMM Regression Results",
          order=c("obstruction",
                  "sdGov",
                  "sdOpp",
                  "govOppDiff",
                  "workload",
                  "change_control",
                  "partyStrength",
                  "maj_size",
                  "daysUntilElect",
                  "reformprevsession",
                  "(Intercept)"),
          covariate.labels = c("obstruction (H1)",
                               "gov. polarisation (H2)",
                               "opp. polarisation (H3)",
                            
                               "gov.-opp. polarisation (H4)",
                               
                               
                               
                               "workload",
                               "party control change",
                               "party strength",
                               "majority",
                               "days until election",
                               "reform t - 1",
                               "intercept"),
          #column.labels=c("coefficient","SE"),
          single.row=TRUE,
          dep.var.caption  = "Dependent variable:",
          dep.var.labels   = "Anti-dilatory reform (x = 1), or not (x = 0)",style="apsr",
          report="vc*s",float=T)

######################
#Split duration models
######################
library("spduration")
UKHCSOdat$state_opening <- as.Date(UKHCSOdat$state_opening)
UKHCSOdat <- add_duration(UKHCSOdat, "dv_bn", unitID = "parliament_ref", tID = "session_number",
                          freq = "year", ongoing = FALSE)


#model 1
weib_model1 <- spdur(
  
  #duration
  duration ~  
    sdGov +
    govOppDiff +
    maj_size + 
    obstruction +
    partyStrength +
    change_control + 
    workload +
    daysUntilElect + 
    reformprevsession,
  
  
  #risk
  atrisk ~ 
    sdGov +
    govOppDiff +
    maj_size + 
    obstruction +
    partyStrength +
    change_control + 
    workload +
    daysUntilElect + 
    reformprevsession,
  data = UKHCSOdat, distr = "weib", silent = TRUE)

#model 2
weib_model2 <- spdur(
  
  #duration
  duration ~  
    sdOpp +
    maj_size + 
    obstruction +
    partyStrength +
    change_control + 
    workload +
    daysUntilElect + 
    reformprevsession,
  
  
  #risk
  atrisk ~ 
    sdOpp +
    maj_size + 
    obstruction +
    partyStrength +
    change_control + 
    workload +
    daysUntilElect + 
    reformprevsession,
  data = UKHCSOdat, distr = "weib", silent = TRUE)

# generate and write out regression table
par(mfrow=c(1,1))
source("customFunctions/genRegTableSpdur.R")

durationResults <- cbind(genRegTableSpdur(weib_model1)$durationResults,
                genRegTableSpdur(weib_model2)$durationResults
)

riskResults <- cbind(genRegTableSpdur(weib_model1)$riskResults,
                         genRegTableSpdur(weib_model2)$riskResults
)

fitStatsN <- rbind(cbind(genRegTableSpdur(weib_model1)$fitStats,
                                  genRegTableSpdur(weib_model2)$fitStats
                                 
),cbind(genRegTableSpdur(weib_model1)$obs,
        genRegTableSpdur(weib_model2)$obs)
)

#write out fit plots
pdf(file="figures/weib_model1SepPlot.pdf",height=3,width=3)
plot(weib_model1)
dev.off()

pdf(file="figures/weib_model2SepPlot.pdf",height=3,width=3)
plot(weib_model2)
dev.off()

output.file <- "models/splitDurModelsDuration.tex"
write("",output.file)
for(i in 1:nrow(durationResults)){
  cat(paste(gsub("NA","",paste(rownames(durationResults[i,]),paste(paste(unlist(durationResults[i,]),collapse = " & "),"\\\\",sep=" "),sep=" & "))),"\n",file = output.file,append=T)

}

output.file <- "models/splitDurModelsRisk.tex"
write("",output.file)
for(i in 1:nrow(riskResults)){
  cat(paste(gsub("NA","",paste(rownames(riskResults[i,]),paste(paste(unlist(riskResults[i,]),collapse = " & "),"\\\\",sep=" "),sep=" & "))),"\n",file = output.file,append=T)
  
}

output.file <- "models/splitDurModelsFitStats.tex"
write("",output.file)
for(i in 1:nrow(fitStatsN)){
  cat(paste(paste(rownames(fitStatsN[i,]),paste(paste(unlist(fitStatsN[i,]),collapse = " & "),"\\\\",sep=" "),sep=" & ")),"\n",file = output.file,append=T)
  
}


#############
#generate coefficient plots
statSigCol <- "black"
statInsignCol <- "grey80"
ratioFormat <- T
xlab <- "Log odds"
source("customFunctions/genCoefplotSpldur.R")
mod1Plot <- genCoefplotSpldur(mod = weib_model1,
                        leftAxis = T,
                        statSigCol = statSigCol,
                        statInsignCol = statInsignCol,
                        plotTitle = "Model 1",
                        xlab = xlab)

mod2Plot <- genCoefplotSpldur(mod = weib_model2,
                              leftAxis = F,
                              statSigCol = statSigCol,
                              statInsignCol = statInsignCol,
                              plotTitle = "Model 2",
                              xlab = xlab)


library(grid)
pdf(file="figures/FIG2_SPLDUR_COEFPLOT1.pdf",width =.8*8, height = .8*3)
grid.newpage()
grid.draw(cbind(ggplotGrob(mod1Plot$durPlot), 
                ggplotGrob(mod2Plot$durPlot), 
                
                size = "last"))
dev.off()


library(grid)
pdf(file="figures/FIG2_SPLDUR_COEFPLOT2.pdf",width =.8*8, height = .8*3.5)
grid.newpage()
grid.draw(cbind(ggplotGrob(mod1Plot$riskPlot), 
                ggplotGrob(mod2Plot$riskPlot), 
                size = "last"))
dev.off()

#############
#effect plots
#############
#high risk
pdf(file="figures/FIG3_LOGTIME_SDGOV.pdf",width =.8*8, height = 3.5)
par(mfrow=c(1,2))
plot_hazard(weib_model1,
            xvals = with(UKHCSOdat,c(1, #sdGov 2 sd below mean                        
                                     sdGov = -1,
                                     mean(govOppDiff,na.rm=T),
                                     
                                     mean(maj_size,na.rm=T),
                                     mean(obstruction,na.rm=T),
                                    mean(partyStrength,na.rm=T),
                                     mean(change_control,na.rm=T),
                                     mean(workload,na.rm=T),
                                     mean(reformprevsession,na.rm=T),
                                     mean(daysUntilElect,na.rm=T)
            )),
            zvals = with(UKHCSOdat,c(1,mean(sdGov,na.rm=T),
                                     mean(govOppDiff,na.rm=T),
                                     mean(maj_size,na.rm=T),
                                     obstruction = 1, #obstruction 2 sd above mean
                                    mean(partyStrength,na.rm=T),
                                     mean(change_control,na.rm=T),
                                     mean(workload,na.rm=T),
                                     mean(reformprevsession,na.rm=T),
                                     mean(daysUntilElect,na.rm=T)
            )),
            ci=T,
            t=c(1:5),
            main = "gov. polarisation (high risk)",
            cex.main=0.75)



#low risk
plot_hazard(weib_model1,
            xvals = with(UKHCSOdat,c(1,
                                     sdGov = 1, #sdGov 2sd above mean
                                      govOppDiff = 1,
                                     mean(maj_size,na.rm=T),
                                     mean(obstruction,na.rm=T),
                                     mean(partyStrength,na.rm=T),
                                     mean(change_control,na.rm=T),
                                     mean(workload,na.rm=T),
                                     mean(reformprevsession,na.rm=T),
                                     mean(daysUntilElect,na.rm=T)
            )),
            zvals = with(UKHCSOdat,c(1,mean(sdGov,na.rm=T),
                                     mean(govOppDiff,na.rm=T),
                                     mean(maj_size,na.rm=T),
                                     obstruction = -1, #obstruction 2sd below mean
                                    mean(partyStrength,na.rm=T),
                                     mean(change_control,na.rm=T),
                                     mean(workload,na.rm=T),
                                     mean(reformprevsession,na.rm=T),
                                     mean(daysUntilElect,na.rm=T)
            )),
            ci=T,
            t=c(1:5),
            main = "gov. polarisation (low risk)",
            cex.main=.75)

dev.off()


#govOppDiff
pdf(file="figures/FIG4_LOGTIME_GOVOPPDIFF.pdf",width =.8*8, height = 3.5)
par(mfrow=c(1,2))

#high risk
plot_hazard(weib_model1,
            xvals = with(UKHCSOdat,c(1,                                     mean(sdGov,na.rm=T),
                                     govOppDiff = -1,  #1sd below mean
                                     
                                     mean(maj_size,na.rm=T),
                                     mean(obstruction,na.rm=T),
                                    mean(partyStrength,na.rm=T),
                                     mean(change_control,na.rm=T),
                                     mean(workload,na.rm=T),
                                     mean(reformprevsession,na.rm=T),
                                     mean(daysUntilElect,na.rm=T)
            )),
            zvals = with(UKHCSOdat,c(1,mean(sdGov,na.rm=T),
                                     mean(govOppDiff,na.rm=T),
                                     mean(maj_size,na.rm=T),
                                     obstruction = 1,
                                     mean(partyStrength,na.rm=T),
                                     mean(change_control,na.rm=T),
                                     mean(workload,na.rm=T),
                                     mean(reformprevsession,na.rm=T),
                                     mean(daysUntilElect,na.rm=T)
            )),
            ci=T,
            t=c(1:5),
            main = "opp. polarisation (high risk)",
            cex.main=.75)



#low risk
plot_hazard(weib_model1,
             xvals = with(UKHCSOdat,c(1,
                                      mean(sdGov,na.rm=T),
                                      govOppDiff = 1,
                                      mean(maj_size,na.rm=T),
                                      mean(obstruction,na.rm=T),
                                      mean(partyStrength,na.rm=T),
                                      mean(change_control,na.rm=T),
                                      mean(workload,na.rm=T),
                                      mean(reformprevsession,na.rm=T),
                                      mean(daysUntilElect,na.rm=T)
             )),
             zvals = with(UKHCSOdat,c(1,mean(sdGov,na.rm=T),
                                      mean(govOppDiff,na.rm=T),
                                      mean(maj_size,na.rm=T),
                                      obstruction = -1,
                                     mean(partyStrength,na.rm=T),
                                      mean(change_control,na.rm=T),
                                      mean(workload,na.rm=T),
                                      mean(reformprevsession,na.rm=T),
                                      mean(daysUntilElect,na.rm=T)
             )),
             ci=T,
             t=c(1:5),
             main = "opp. polarisation (low risk)",
             cex.main=.75)

dev.off()

#sdOpp
#govOppDiff
pdf(file="figures/FIG5_LOGTIME_SDOPP.pdf",width =.8*8, height = 3.5)
par(mfrow=c(1,2))

#high risk
plot_hazard(weib_model2,
            xvals = with(UKHCSOdat,c(1,
                                     sdOpp = -1,
                                   
                                     mean(maj_size,na.rm=T),
                                     mean(obstruction,na.rm=T),
                                     mean(partyStrength,na.rm=T),
                                     mean(change_control,na.rm=T),
                                     mean(workload,na.rm=T),
                                     mean(reformprevsession,na.rm=T),
                                     mean(daysUntilElect,na.rm=T)
            )),
            zvals = with(UKHCSOdat,c(1,mean(sdOpp,na.rm=T),
                                     
                                     mean(maj_size,na.rm=T),
                                      obstruction = 1,
                                     mean(partyStrength,na.rm=T),
                                     mean(change_control,na.rm=T),
                                     mean(workload,na.rm=T),
                                     mean(reformprevsession,na.rm=T),
                                     mean(daysUntilElect,na.rm=T)
            )),
            ci=T,
            t=c(1:5),
            main = "gov-opp. polarisation (high risk)",
            cex.main=.75)



#low risk
plot_hazard(weib_model2,
             xvals = with(UKHCSOdat,c(1,
                                      sdOpp = 1,
                                      
                                      mean(maj_size,na.rm=T),
                                      mean(obstruction,na.rm=T),
                                      mean(partyStrength,na.rm=T),
                                      mean(change_control,na.rm=T),
                                      mean(workload,na.rm=T),
                                      mean(reformprevsession,na.rm=T),
                                      mean(daysUntilElect,na.rm=T)
             )),
             zvals = with(UKHCSOdat,c(1,mean(sdGov,na.rm=T),
                                  
                                      mean(maj_size,na.rm=T),
                                      obstruction = -1,
                                      mean(partyStrength,na.rm=T),
                                      mean(change_control,na.rm=T),
                                      mean(workload,na.rm=T),
                                      mean(reformprevsession,na.rm=T),
                                      mean(daysUntilElect,na.rm=T)
             )),
             ci=T,
             t=c(1:5),
             main = "gov-opp. polarisation (low risk)",
             cex.main=.75)

dev.off()

###################
# Outlier analysis
###################
## Mixed effects model

#Step 1: hold out entire parliaments
# SPLDUR MODELS
#mod 1
source("customFunctions/parlOutlierAnalysis.R")
model <- weib_model1
varName <- "sdGov"

problematicParls <- parlOutlierAnalysis(model,
                                        "SPDUR",
                                        varName)

nrow(problematicParls)
min(problematicParls$pVal)
max(problematicParls$pVal)

#mod 1
model <- weib_model1
varName <- "govOppDiff"

problematicParls <- parlOutlierAnalysis(model,
                                        "SPDUR",
                                        varName)

nrow(problematicParls)
min(problematicParls$pVal)
max(problematicParls$pVal)

#mod 2
model <- weib_model2
varName <- "sdOpp"

problematicParls <- parlOutlierAnalysis(model,
                                        "SPDUR",
                                        varName)

nrow(problematicParls)
min(problematicParls$pVal)
max(problematicParls$pVal)

################
# PER SESSION
################
source("customFunctions/sessionOutlierAnalysis.R")

# SPLDUR MODELS
#mod 1
model <- weib_model1
varName <- "sdGov"

problematicSessions <- sessionOutlierAnalysis(model,
                                              "SPDUR",
                                              varName)



nrow(problematicSessions)
min(problematicSessions$pVal)
max(problematicSessions$pVal)

#mod 1
model <- weib_model1
varName <- "govOppDiff"

problematicSessions <- sessionOutlierAnalysis(model,
                                              "SPDUR",
                                              varName)

nrow(problematicSessions)
min(problematicSessions$pVal)
max(problematicSessions$pVal)

#mod 2
model <- weib_model2
varName <- "sdOpp"

problematicSessions <- sessionOutlierAnalysis(model,
                                              "SPDUR",
                                              varName)

nrow(problematicSessions)
min(problematicSessions$pVal)
max(problematicSessions$pVal)

####################
#Breakpoint analysis
####################

source("customFunctions/breakPointAnalysis.R")
breakPointAnalysis("obstruction",UKHCSOdat)

